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Abstract 

We  propose  stabilized  interior  penalty  discontinuous  Galerkin  methods  for  the  indefinite 
time-harmonic  Maxwell  system.  The  methods  arc  based  on  a  mixed  formulation  of  the 
boundary  value  problem  chosen  to  provide  control  on  the  divergence  of  the  electric  field. 
We  prove  optimal  error  estimates  for  the  methods  in  the  special  case  of  smooth  coefficients 
and  perfectly  conducting  boundary  using  a  duality  approach. 
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1  Introduction 


The  numerical  solution  of  the  time-harmonic  Maxwell’s  equations  presents  a  num¬ 
ber  of  challenges.  First,  away  from  boundaries  and  material  interfaces,  the  solution 
is  smooth  and  oscillatory.  The  need  to  approximate  the  oscillations  requires  a  suffi¬ 
ciently  fine  grid  compared  to  the  wave-length  of  the  solution,  and  results  in  a  large 
number  of  degrees  of  freedom  if  many  wavelengths  are  contained  in  the  domain  of 
interest.  This  requirement  can  be  loosened  (but  not  entirely  avoided)  by  the  use  of 
high  order  methods  [1,2],  so  that  it  is  desirable  to  use  high  order  methods  where 
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the  solution  is  smooth.  A  second  problem  is  that  at  the  boundary  of  the  domain 
the  solution  can  be  singular  [3].  Indeed,  on  a  non-convex  polyhedral  domain  the 
straightforward  application  of  continuous  finite  element  methods  can  result  in  a  dis¬ 
crete  solution  that  converges  to  a  vector  function  that  is  not  a  solution  of  Maxwell’s 
equations  [4].  While  it  is  possible  to  modify  the  variational  form  to  correct  for  this 
failure  [5],  a  similar  problem  also  occurs  at  interfaces  between  different  materials. 
This  is  because  at  discontinuities  in  the  electric  properties  of  the  materials  in  the 
domain  of  the  electromagnetic  field,  the  electric  field  is  discontinuous.  Thus  con¬ 
tinuous  elements  need  to  be  modified  at  such  interfaces.  A  third  problem,  which 
we  will  not  discuss  here,  is  the  numerical  solution  of  the  sparse  indefinite  matrix 
problem  resulting  from  the  finite  element  discretization. 

Considerations  of  the  first  two  problems  mentioned  above  have  lead  to  a  wide¬ 
spread  adoption  of  edge  finite  elements  [6,7]  for  the  discretization  of  the  time- 
harmonic  Maxwell’s  equations.  For  an  engineering  view  of  such  elements,  a  good 
summary  is  contained  in  the  books  [8,9].  An  error  analysis  of  these  elements  has 
been  given  in  [10-12]  and  the  profound  connection  between  these  elements  and 
differential  forms  has  been  noted  for  example  in  [13-15].  Perhaps  the  main  problem 
with  such  elements  is  that  they  become  rather  complex  as  the  order  of  the  elements 
is  increased,  and  like  all  conforming  methods  they  require  a  suitable  finite  element 
grid  which  complicates  implementing  adaptive  solvers.  Nevertheless,  adaptive  hp- 
finite  element  solvers  have  been  implemented  and  show  considerable  promise  [16]. 

In  this  paper  we  propose  a  new  way  to  discretize  the  indefinite  time-harmonic 
Maxwell  system  based  on  a  discontinuous  Galerkin  method  (denoted  DG  in  the 
remainder  of  the  paper).  In  particular,  we  propose  a  suitable  extension  of  the  in¬ 
terior  penalty  methods  to  the  Maxwell  system.  These  methods  date  back  at  least 
to  [17-20]  and  have  been  studied  for  coercive  elliptic  and  convection  diffusion 
problems  more  recently  in  [21-23].  We  mention  that  several  other  DG  methods  for 
standard  coercive  elliptic  problems  can  be  found  in  the  literature  (for  instance  the 
LDG  method  [24,25]  or  the  DG  method  introduced  by  Baumann  and  Oden  [26,27]) 
and  unified  analyzes  of  discontinuous  methods  in  the  context  of  elliptic  problems 
have  been  presented  in  [28,29].  For  the  time-harmonic  Maxwell  equations  in  the 
low-frequency  regime,  where  the  resulting  bilinear  forms  are  coercive,  the  LDG 
methods  have  been  recently  investigated  in  [30]. 

Our  goal  is  to  produce  a  flexible  solver  in  which  the  order  of  the  scheme  can  be 
changed  easily  between  different  regions  of  the  grid.  In  addition  we  hope  to  exploit 
the  fact  that  DG  grids  do  not  need  to  be  aligned  in  order  to  improve  the  efficiency 
of  wave  propagation  of  the  method.  Thus  in  regions  with  different  electromagnetic 
properties  (and  hence  different  wave  speeds),  different  grid  sizes  can  be  used  to  bal¬ 
ance  the  propagation  accuracy  of  the  scheme  in  each  subdomain.  Finally  it  may  be 
possible  to  “tune”  parameters  in  the  DG  scheme  to  improve  propagation  accuracy 
(this  is  certainly  possible  in  ID  one  space  dimension!). 
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In  this  paper  we  prove  basic  error  estimates  for  our  proposed  schemes  under  the 
assumption  of  smoothly  varying  material  properties.  This  assumption  is  needed  for 
certain  a-priori  estimates  used  in  the  analysis.  Ultimately  we  hope  to  extend  these 
results  to  more  general  coefficients  and  boundary  conditions,  and  we  detail  the 
formulation  in  these  cases. 

Perhaps  the  closest  approach  in  the  literature  to  the  DG  methods  we  propose  is  the 
ultra  weak  variational  method  of  Cessenat  [31].  While  successful  in  practice,  this 
method  is  still  incompletely  understood  on  a  theoretical  level.  For  example  con¬ 
vergence  is  not  proved  for  the  standard  perfectly  electrically  conducting  boundary 
condition  (or  near  a  singularity),  or  in  general  throughout  the  domain  of  compu¬ 
tation.  However  the  successful  use  of  this  method  is  one  motivation  for  propos¬ 
ing  the  methods  in  this  paper  which  are  convergent  globally  even  in  the  presence 
of  boundary  singularities.  Another  similar  approach  is  the  mortar  finite  element 
method  applied  to  the  Maxwell  equations  [32].  To  our  knowledge,  convergence  has 
not  been  proved  for  this  method  in  the  case  of  wave  propagation.  However  the  suc¬ 
cess  of  this  method  applied  to  low  frequency  eddy  current  problems  (in  which  case 
the  resulting  bilinear  forms  are  coercive)  suggests  that  mortar  methods  or  similar 
domain  decomposition  methods  could  be  useful  for  scattering  problems  [33].  An¬ 
other  domain  decomposition  approach  is  the  FETI  method  applied  to  the  Maxwell 
equations  [34].  A  Lagrange  multiplier  based  version  of  this  method  was  analyzed 
in  [35]  for  the  coercive  Maxwell  problem  arising  in  time  stepping.  Again  to  our 
knowledge,  convergence  has  not  been  proved  for  this  method  in  the  case  of  wave 
propagation  in  Maxwell’s  equations. 

The  outline  of  our  paper  is  as  follows.  In  section  2,  we  start  by  detailing  the  mixed 
formulation  we  shall  use  as  the  basis  of  the  DG  methods  proposed  here.  We  also 
summarize  some  regularity  and  existence  results.  Then  in  section  3  we  propose 
the  DG  methods  that  are  the  subject  of  this  paper.  The  main  result  of  the  paper 
is  an  optimal  a  priori  error  bound  that  we  present  in  section  4.  Its  proof  is  based 
on  a  duality  approach  and  is  contained  in  section  5  and  section  6.  We  end  our 
presentation  with  some  concluding  remarks  in  section  7. 


2  A  mixed  formulation  for  the  time-harmonic  Maxwell  equations 


In  this  section,  we  introduce  the  time-harmonic  Maxwell  equations  and  present  a 
mixed  formulation  for  the  continuous  problem  which  will  be  the  basis  for  the  DG 
methods  introduced  here. 
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2. 1  Time-harmonic  Maxwell ’s  equations 


Let  be  a  bounded  Lipschitz  polyhedron  in  R3  with  connected  boundary  <9fL  The 
model  problem  we  shall  consider  is  to  compute  a  time-harmonic  electric  field  £ 
in  the  cavity  Q  with  perfectly  conducting  boundary.  Let  co  denote  the  temporal 
frequency  of  the  time-harmonic  field  so  that  the  corresponding  time  dependent  field 
E  at  position  xeSl  and  time  t  is  given  by 

E(x,  t)  =  5ft  (£(x)  exp(— iut)) . 

Then  £  :  — >■  C3  satisfies  the  Maxwell  system 

V  x  jUyT1  V  x  £  —  k2er£  =  J  in  Q,  (1) 

where  //T  is  the  relative  magnetic  permeability  and  er  is  the  relative  electric  permit¬ 
tivity  of  the  medium  in  the  cavity  Q.  We  assume  that  /<T  and  er  are  real,  smooth,  and 
uniformly  positive  functions  of  the  position  in  Q.  In  addition  the  real  wave  number 
k  is  given  by 

k  =  U^/£0fi0, 

where  /i0  is  the  magnetic  permeability  and  e0  is  the  electric  permittivity  of  free 
space.  The  source  function  J  is  related  to  the  applied  current  density  driving  the 
cavity  and  is  assumed  to  be  a  given  vector  function  in  L2( f))3. 

The  assumption  that  Q  has  a  perfectly  conducting  boundary  gives  the  following 
boundary  condition  on  <90: 


n  x  £  =  0  on  <9f I.  (2) 

Here  n  denotes  the  outward  normal  unit  vector  to  <90. 

Throughout  the  paper,  we  will  assume  that  k2  is  not  an  interior  Maxwell  eigenvalue 
(see  also  Proposition  1  below),  i.e.,  for  any  £  ^  0,  the  pair  (A  =  k2,  £)  is  not  an 
eigensolution  of  the  problem  V  x  x  £  =  \er£  in  O,  n  x  £  =  0  on  <90. 

Note  that  this  assumption  would  not  be  necessary  if  some  region  of  O  (containing 
a  ball  of  non-zero  radius)  had  a  non  zero  conductivity  which  would  imply  that  the 
imaginary  part  of  er  is  positive  there.  Note  also  that  in  the  special  case  considered 
here  the  real  and  imaginary  parts  of  the  solution  decouple,  and  hence  we  can  assume 
that  £  is  real.  If  er  is  complex  valued  or  if  impedance  boundary  conditions  are 
imposed,  the  real  and  imaginary  parts  are  coupled. 


2.2  Mixed  formulation 


Our  DG  method  is  based  on  a  mixed  formulation  of  the  Maxwell  boundary  value 
problem  (l)-(2).  Such  formulations  have  been  used  previously  for  edge  element 
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discretizations  of  Maxwell’s  equations  to  improve  stability  [16],  and  to  handle  co¬ 
ercive  problems  in  which  meshes  are  not  aligned  at  a  material  boundary  [35].  We 
can  derive  this  formulation  by  using  a  Helmholtz  decomposition. 

Given  a  domain  D  in  R2  or  R3,  we  denote  by  Hs(D)d,  d  =  1,2,  3,  the  Sobolev 
space  of  real  or  complex  scalar-  or  vector- valued  functions  with  regularity  exponent 
s  >  0,  endowed  with  the  usual  norm  ||  •  \\SjD  and  seminorm  |  •  \SjD.  We  write  Hq(D) 
for  the  subspace  of  Hl(D)  of  functions  with  zero  trace  and  set  L2(D)d  =  H°(D)d. 
For  the  computational  domain  Q  C  R3 ,  we  let 

H( div°r;H)  =  ju  £  L2(H)3  |  V  •  (eyu)  =  0  in  flj 
H( div°;  Q)  =  ju  £  L2(H)3  |  V  •  u  =  0  in  , 

equipped  with  the  L2  (£2) 3 -norm.  We  shall  also  need  to  use  the  standard  spaces 

H(  curl;H)  =  ju  e  L2(H)3  VxuG  L2(H)3| 

H0(  curl;  O)  =  {u  G  if  (curl;  Q)  n  x  u  =  0  on  30}  , 

endowed  with  the  norm 


IMlLi,n  =  llullo,o  +  llv  x  ullo,m 

Furthermore,  let  L2r{ H)3  denote  the  space  of  square  integrable  functions  on  Q 
equipped  with  the  inner  product 

(u,  v)£  =  /  £rU  ■  vdx. 

Jfl 

We  consider  the  Ll2r  (H)3-orthogonal  Helmholtz  decomposition  of  the  function  £  £ 
H0( curl;  H)  given  by  £  =  u  +  Vp,  where  u  £  iF0(curl;  H)  n  H( div°r;  H)  and 
p  £  Hq(Q),  see  [36]. 

Using  this  decomposition,  problem  (l)-(2)  can  be  reformulated  as  follows:  find 
u  £  H (curl;  H)  and  p  £  fF1(H)  such  that 


V  x  /i^V  xu  —  k2eru  —  k2£r'Vp  =  J  in  Q  (3) 

V  •  (eru)  =  0  in  H  (4) 

n  x  u  =  0  on  d£l  (5) 

p  =  0  on  dQ.  (6) 


We  start  by  showing  well-posedness  of  problem  (3)-(6). 

Proposition  1  Assume  that  k 2  is  not  a  Maxwell  eigenvalue.  Then  problem  (3 )-( 6) 
has  a  unique  solution  (u, p)  £  iF0(curl;  H  iF(div°r;  H)  x  Hq(Q),  with  /q^V  x 
u  £  H (curl;  Q),  and  we  have  the  stability  estimates 

||u||curl;n  ^  Cstabll^llo,^;  ||p||l,n  ^  k  Cell ||  3 1 1 0,0 5 
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with  positive  constants  Cstab  cmd.  Ceii,  Ce n  independent  of  k. 


Proof  Consider  the  Helmholtz  decomposition  of  J  as  a  function  in  L2(f))3,  J  = 
Jo  +  Vj  with  J0  G  ff(div°;  fi)  and  j  G  ffg(H).  Owing  to  the  orthogonality  of 
this  decomposition,  problem  (3)-(6)  decouples  into  two  independent  subproblems, 
namely  into  the  Maxwell  problem  with  divergence  free  data 

V  x  pf1  V  xu—  k2sr  u  =  J0  in  Q 

V  •  (sru)  =  0  in  0  (7) 

n  x  u  =  0  on  dQ, 

and  the  elliptic  problem 

—k2V  •  (£rVp)  =  F  in  O,  p  =  0  on  90,  (8) 

with  right  hand  side  F  G  77_1(0)  defined  by  F(q)  =  —  Jn  Vj  ■  Vq  dx,  for  all 
q  G  Hq(Q).  Existence  and  uniqueness  of  solutions  to  (7)  follow  now  in  a  standard 
way  from  Fredholm  theory  and  the  coercivity  of  the  form  J^pf1V  x  u  •  V  x  v  dx 
on  the  space  H0( curl;  0)  n  ff(div°  ;  0)  (see  [36,  Proposition  2.7]).  Furthermore, 
pfxV  x  u  G  if  (curl;  0)  and 


||u||Curl,n  5:  C*stab ||  Jo  ||o,n  5:  C^stab 1 1  1 1 

with  a  stability  constant  Cstab  >  0  depending  on  0,  pr,  er  and  on  the  wave  number 
k2.  For  problem  (8),  existence  and  uniqueness  follows  from  standard  elliptic  theory, 
and  we  have 

IMkn  <  /c-2Ceii||F||_l!n  <  Ar2Cell||V.?||o,n  <  Ar2C'eu||.7||o,n, 

with  Ceii  >  0  only  depending  on  Q  and  er .  □ 

In  our  duality  approach  in  section  6  we  will  also  make  use  of  the  following  regu¬ 
larity  result. 

Proposition  2  For  smooth  coefficients  pr  and  er,  there  exists  a  regularity  exponent 
o  =  cr(fi)  >  |  such  that  the  solution  u  in  (3 )-( 6)  satisfies  u  G  Ha( fl)3  and 
V  x  u  G  Ha(Q)3.  Furthermore, 

||u||CTjn  +  ||  V  x  u||a!n  Creg||^7"||o)o; 

with  a  positive  constant  Creg  depending  on  If2,  pr  and  er. 

Proof  By  decoupling  problem  (3)-(6)  into  (7)  and  (8),  we  see  that  the  solution  u 
of  (7)  satisfies  V  x  u  G  L2(fi)3,  V  •  u  G  L2( fi)3  (here  we  use  the  assumption  that 
er  is  smooth,  so  that  erV  ■  u  =  V  •  (eTu)  —  Ver  •  u  holds  true)  and  n  x  u  =  0 
on  dfl  .  From  [37,  Proposition  3.7],  it  follows  that  u  G  Hai(fl)3  for  <7i  >  \  and 
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||u||fflin  <  C'(Q)||u||curi)n  with  an  embedding  constant  C(Q)  just  depending  on  fi. 
Thus,  from  Proposition  1,  we  have  HuH^n  <  C1!  ||  J7'||0,o- 

Now  set  w  =  x  u.  From  the  first  equation  in  (7),  we  have  V  x  w  =  J0  + 

k2er  u  G  L2(Q)3.  Furthermore,  V  •  (pr  w)  =  0  and  pr  w  -  n  =  Vxu-n  =  0on  dQ. 
Since  pr  is  smooth,  using  again  [37,  Proposition  3.7],  it  follows  that  w  G  Ha2(Q.)3 
for  a2  >  |  and  ||w||ff2in  <  C^jllwHcuri^  <  C(fl,  pr){\\ V  x  u||g  n  +  ||J0  + 
A;2£ru||o:n)  2 .  Hence,  from  the  triangle  inequality  and  Proposition  1,  we  conclude 
that  ||w||a2!n  <  C2II  Jl|o,n-  Choosing  a  :=  min{cri,cr2}  and  Creg  =  rriax{C'iX'2} 
completes  the  proof.  □ 

Remark  3  If  the  polyhedron  Q  is  convex  and  pr  =  er  =  1,  the  parameter  o  in 
Proposition  2  can  be  chosen  as  o  =  1,  see  [37]. 


3  Discontinuous  Galerkin  discretization 


In  this  section,  we  introduce  stabilized  interior  penalty  discontinuous  Galerkin  dis¬ 
cretizations  for  the  Maxwell  system  (3)-(6). 


3. 1  Triangulations 


Let  Th  be  a  regular  triangulation  of  the  domain  Q  into  tetrahedra.  We  denote  by  hK 
the  diameter  of  the  element  K  and  set  h  =  m&xKeph  hK.  The  diameter  of  the  face 
/  is  denoted  by  hj.  We  also  assume  the  triangulation  to  be  shape  regular,  that  is, 
there  is  a  positive  constant  k  such  that,  for  any  K  G  Th, 

hx  . 

—  <  K, 

Pk 

where  pK  is  the  diameter  of  the  biggest  ball  contained  in  K  (see  [38,  p.  124]). 

Let  S  be  the  union  of  all  the  faces  of  Th,  and  Ex  the  union  of  the  internal  faces.  We 
define  the  function  h  in  L°°(E )  by 


h  =  h(x)  =  hf  if  x  G  /. 


3.2  Trace  operators 


First,  we  need  to  define  some  notation  concerning  functions  in  Hs{Th )  :=  {v  : 
v\K  g  Hs(K),  I<  g  Th},  for  s  >  The  elementwise  traces  of  such  functions 
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belong  to  TR(£)  :=  UKerhL2(dK)-,  they  are  double-valued  on  Ex  and  single¬ 
valued  on  E  \  St-  The  space  L2(E )  can  be  identified  with  the  functions  in  TR(£) 
for  which  the  two  trace  values  coincide. 

Next,  we  introduce  certain  trace  operators.  To  this  end,  fix  w  G  TR(£)3  and  <p  G 
TR(£),  and  let  e  C  Ex  be  an  interior  face  shared  by  the  elements  K\  and  I<2-  Let 
rq  be  the  normal  unit  vector  pointing  exterior  to  Kt  and  w,  =  w\gKi,  (pi  =  p\dKt 
(i  =  1.  2).  Then  we  define  for  x  G  e  the  average,  the  tangential  jump  and  the 
normal  jump  of  w  as  follows: 

]{w}}  =  ^(wi+w2)  [w]T  =  mxw1+n2xw2  [w]N  =  wi-m+w2-n2. 
Similarly,  we  define  for  x  G  e  the  average  and  the  normal  jump  of  p  by 

=  2^  +  ^2)  Miv  =  <Pini  +  p2n2- 
Then,  on  any  boundary  face  e  C  S  \  Ex-,  we  set  for  x  G  e 

Hwlf  =  w  Mr  =  n  x  w  Miv  =  pn. 

Since  we  will  not  require  either  of  the  quantities  j{ and  |w]  N  on  the  boundary 
E  \  Ex,  we  leave  them  undefined. 

If  w  G  H (curl;  Q),  then,  for  all  e  C  Ex,  the  jump  condition  ni  x  wi  +  n2  x  w2  =  0 

holds  true  in  H00 2  (e)3,  and  thus  also  in  L2(e)3  (for  the  definition  of  H00 2  (e),  see, 
e.g.,  [39]).  Therefore,  [w]T  is  equal  to  zero  on  Ex.  Similarly,  for  w  e  H( div;  fl), 
we  have  that  [w]jv  is  well-defined  and  equal  to  zero  on  Ex-  Furthermore,  for  the 
exact  solution  u  G  ff0(curl;  El)  n  H( div°r ;  fl),  we  have  [u]T  =  0  in  L2(e)3  for  any 
boundary  face  e,  in  addition  to  [u]T  =  0  and  [cru];V  =  0  on  £x. 


3.3  Stabilized  interior  penalty  discontinuous  Galerkin  methods 


We  approximate  u  and  p  in  the  discontinuous  finite  element  space  V/t  x  Qh  where 

V/i  =  jv  G  L2(Elf  :  v|*  G  \/K  G  Th} 

Qh  =  {qe  L2(Q)  :  q\K  G  V\K\  VK  G  %}, 

for  an  approximation  order  t  >  1,  with  Tl(K)  denoting  the  space  of  polynomials 
of  degree  at  most  l  on  K. 

We  consider  the  following  discontinuous  Galerkin  method:  find  (u/i,  ph )  G  SI  hxQh 
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such  that,  for  any  (v,  q)  e  \h  x  Qh, 


(9) 

(10) 


a(uh,v)+c(uh,v)+d(uh,v)-k2(uh,v)Er+b(v,ph)  =  :F(v) 

b(uh,q)  ~e(ph,q)  =0. 


Here, 


a(u,  v) 

b(v,p) 
c(u,  v) 
d(u.  v) 


e(p,  q) 


x  u-  Vft  x  v  dx  J  [u]t  ■  "g/A-  1 V h  x  v}}  ds 


~  f£  Iv1t  •  lpT  1 V.  x  ufds 

=/c2  f  p  S7h  •  (erv)  dx  -  k2  f  -gp}}  |erv]jV  ds 
Jo 

=a;  J  h_1m_1[u]r  •  [v]T  ds 

=/c2/3  /  h  [£ru]iv  [erv]jv  ds 
JSx 

+  A;2/?  hK  [  Vft  ■  (eru)  VA  •  (erv)  dx 
=fc27jf  h_1e[p]jv  •  [§]Ards  +  a  J  h_1m_1[Vftp]r  •  [Vftg]Tds, 


with  V/jX,  V/j-  and  V/t  denoting  the  elementwise  curl,  divergence  and  gradient, 
respectively,  and  the  functions  m  and  e  are  defined  on  £  as  the  restriction  to  S  of 
pr  and  er,  respectively.  The  parameters  a,  [3  and  7  in  the  forms  c,  d,  and  e  are 
positive.  The  purpose  of  these  forms  is  to  stabilize  the  method.  The  functional  T  at 
right-hand  side  of  (9)  is 

3F{y)  =  [  J  •  v  dx. 

Jn 

Let  us  discuss  the  following  points  about  this  method: 


•  The  form  a(-,  ■)  +  c(-,  •)  corresponds  to  the  interior  penalty  discretization  of  the 
curl-curl  operator,  cf.  [28];  it  is  symmetric  and  stable  provided  that  the  parameter 
a  is  large  enough  (see  Lemma  14  below).  The  nonsymmetric  variant  of  the  interior 
penalty  discretization  is  obtained  by  replacing  a  by 


(u,  v)  =|/ir1ViXU'Vtx  vdx  +  |  [u]T  •  {{pr  x  v} 

-  J£Mt  ■  iPr^h  X  U }}ds. 


ds 


(ID 


Then  the  form  a(-,  •)  +  c(-,  •)  is  nonsymmetric,  but  stable  for  any  a  >  0  (see  Re¬ 
mark  15  below).  In  the  following  we  will  only  present  the  analysis  for  the  symmet¬ 
ric  method  in  (9)-(10),  but  emphasize  that  the  error  estimates  so  obtained  hold  true 
verbatim  for  its  nonsymmetric  variant. 


•  The  form  &(•,•)  discretizes  the  divergence  constraint  in  the  mixed  formulation 
(3)-(6)  by  means  of  DG  techniques,  similar  to  the  forms  used  in  [40]  for  the  Stokes 
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system.  Notice  that,  after  integration  by  parts,  the  form  b(-,  •)  can  also  be  expressed 
by 

&(v,p)  =  -k2  ■  Vhp  dx  +  k2  J  \p]N  •  (v,p)  G  Vft  x  Qh. 


•  The  forms  df.  •)  and  e(-,  •)  provide  stabilization.  While  df,  •)  is  related  to  the 
divergence  constraint,  the  form  e(-,  •)  provides  stability  via  control  of  jumps  of  the 
scalar  potential  p.  We  found  it  necessary  to  include  these  forms  in  order  to  be  able 
to  prove  optimal  error  estimates  with  our  techniques  of  analysis.  Whether  or  not 
similar  results  can  actually  be  obtained  without  these  stabilization  forms  remains 
an  open  question  and  will  be  investigated  numerically  in  a  forthcoming  work. 


•  By  elementary  manipulations  the  third  term  in  the  form  a(-,  •)  can  be  expressed 
by 


-  J  Wt  •  {{dr  1 X  u}}  ds 


v  •  nK  x  [pr  1  V/t  x  u]  ds 


-  f  W  •  \pr  1 V,  x  u]Tds, 
J  £x 


(12) 


for  all  u,  v  G  V/(.  where  n K  is  the  outward  normal  unit  vector  to  OK.  For  the  exact 
solution  u  we  have  u  G  ff0(curl;  f2)  and  x  u  G  ff(curl;  f2).  Thus,  [u]t  =  0 
on  8  and  fl/t"1  V&  x  u]t  =  0  on  Ex,  and,  with  (12),  o(u,  v)  has  to  be  understood  as 


o(u,  v)  =  /  /i  Vft  XU  •  V/,  x  vdx+  V  /  V  •  nK  x  [p  1  Wh  x  u]  ds , 
J  n  ,77,,  JdK 


K£K  ' 


for  v  G  V/i,  where  the  boundary  integrals  are  in  fact  duality  pairings. 


Let  us  now  address  the  consistency  of  the  method. 


Proposition  4  The  discontinuous  Galerkin  method  in  (9)-(10)  is  consistent,  i.  e. , 
the  exact  solution  (u.  p)  of  problem  (3)-(6)  satisfies  (9)-(10),  for  all  test  functions 
(v,q)  E\hx  Qh. 


Proof  We  have  u  G  ff0(curl;  Q),  //“1V  x  u  G  if  (curl;  fi)  and  u  G  H{  div°r;  fi). 
Thus,  [u]r  =  0  on  8,  as  well  as  [£ru]jv  =  0  and  |/rrT1V/l  x  ujT  =  0  on  8x- 
Moreover,  p  G  (91)  and  Vp  G  H0( curl;  fi)  and  thus  |p]iV  =  0  and  [VpJ'i'  =  0 
on  8.  The  second  equation  (10)  is  then  trivially  satisfied  for  all  q  G  Qh-  From  (12), 
the  first  equation  (9)  reduces  to 


/  /i/Vxu'VftXvrfxf  V  /  v  •  n  K  x  \/j  1'Vh  x  ul  ds 
Jn  JdK 

—  k2  eru  ■  v  dx  +  /c2  /  p  V/j  •  (£rv)  dx.  —  k2  /  jfpll^T'Vjjvds  =  j  J  ■  v  dx, 
Jn  Jn  J£x  Jn 
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for  v  e  V/,,.  Integration  by  parts  over  each  element  K.  taking  into  account  the 
boundary  condition  for  u  and  p,  yields 


V  x  pr  1 V  x  u  —  k2er u  — 


v  dx, 


which  is  satisfied  for  all  v  e  V/j. 


□ 


Remark  5  Note  that  the  regularity  of  the  solution  stated  in  Proposition  2  as  con¬ 
sequence  of  the  smoothness  assumption  on  the  coefficients  is  needed  neither  in  the 
definition  of  the  method,  nor  in  the  proof  of  Proposition  4.  As  a  matter  of  fact,  the 
method  is  defined  and  consistent  for  piecewise  smooth  coefficients  pr  and  er.  In 
this  case,  the  functions  m  and  e  have  to  be  adjusted  by  taking  the  corresponding 
averages,  i.  e. ,  by  taking  m  =  and  e  =  {{cr  }}  on  S. 


Remark  6  The  analysis  developed  in  the  following  sections  makes  use  of  con¬ 
forming  projection  operators,  and  therefore  only  covers  the  case  of  meshes  that 
do  not  contain  hanging  nodes.  On  the  other  hand,  the  DG  method  is  well-defined 
for  general  non-matching  grids.  In  this  case,  the  interior  faces  are  understood  as 
the  ( non-empty )  interiors  of  the  intersections  between  two  adjacent  elements  and 
the  function  h  has  to  be  redefined  on  £  as 


h  =  h(x) 


h(hK  +  hKi )  ijx  £  e  =  dK  n  dK' 

< 

hK  ifs.  (E  e  C  dK  n  dfl. 


4  The  main  result 


In  this  section,  we  present  and  discuss  our  main  result  -  an  optimal  a  priori  error 
estimate  for  the  DG  method  in  (9)-(10).  The  proof  of  this  bound  is  developed  in 
section  5  and  section  6;  it  is  based  on  a  suitable  duality  argument  that  heavily  relies 
on  the  regularity  result  of  Proposition  2,  and  therefore  the  assumption  of  smooth 
coefficients  pr  and  eT  is  essential.  Moreover,  in  order  to  simplify  the  presentation, 
we  restrict  ourselves  to  the  case  pr  =  er  =  1.  The  extension  to  general  smooth 
coefficients  is  straightforward. 

Define  the  spaces 

V(h)  :=  Vh  +  (tf0(curl;  Q)  n  H{ div°;  Q))  Q(h)  :=  Qh  + 
and  the  broken  norm  |||  (u,  p)  |||/t  given  by 

III  (u,p)  I*  =  ^2||u||oin  +  *2||Vftp||oin  +  I  (u,p)  1 1, 
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where  the  seminorm  |  (u,  p)  7  is  given  by 

I  (u,p)  1 1  =  ||Vfc  X  u 1 1 o,r2  +  a  ||n  ^Mr||o,f  +  Ilh^Julivllo^ 

+  k2p  ^llV-ullo,iv  +  /c27||h_^b]iv||o)£+a||^^[V/lp]T||o,£- 

KeTh 

It  is  easy  to  see  that  |||  (•,  •)  is  actually  a  norm  in  V(h)  x  Q(h). 

Our  main  result  establishes  error  estimates  in  the  norm  |||  (u,  p )  |||/t. 

Theorem  7  Assume  that  the  exact  solution  (u.  p)  of  the  continuous  problem  (3)- 
(6)  satisfies 

U  e  HS(Q )3  V  X  U  e  HS(Q)3  P  e  HS+1(Q)  s  >  \,  (13) 

£ 

and  let  (u h,Ph)  be  the  discrete  solution  obtained  by  the  DG  method  (9)-(10). 
Then,  there  exists  positive  constants  ao  and  /30,  with  o0  =  afin,  l)  and  /30  = 
/3o(k,  t-  Ceii),  such  that  for  a  >  a0  and  j3  >  (30  we  have  the  error  bound 

I  (u-  u h,p-ph)  ||ft  <  C/imm{As}(||u||S;n  +  ||  V  x  u||a,n  +  ||p||a+i,n), 

provided  that  h  <  h0for  some  h0  =  h0(K,  £,  k,  Creg,  a,  ft,  7,  0,  s). 

Remark  8  Theorem  7  guarantees  optimal  a  priori  error  bounds  provided  that  the 
stabilization  parameters  a  and  (3  are  large  enough.  Restrictions  of  this  type  are 
typically  encountered  in  interior  penalty  methods  (see,  e.g.,  [28,41]).  It  is  worth 
noting  that  ao  and  fio  are  independent  of  the  wave  number  k2. 

Remark  9  Note  that  the  smoothness  assumptions  u  6  HS(Q)  ’  and  V  x  u  6 
Hs  (Q)3  follow  from  Proposition  2,  whereas  the  assumption  p  £  HS+1(Q)  does  not 
seem  to  hold  true  for  general  source  terms  ff  in  L2(Q)3.  This  lack  of  smoothness  in 
the  potential  pfor  general  right  hand  sides  will  be  the  major  difficulty  in  our  duality 
argument  in  section  6  below  and  is  the  reason  why  we  introduced  the  stabilization 
forms  df,  •)  and  e(-,  •).  We  also  point  out  that  for  divergence  free  source  terms  J 
(often  encountered  in  practice )  we  have  p  =  0  and  the  assumption  p  £  HS+1(Q)  is 
trivially  satisfied. 

Proceeding  along  the  lines  of  [42],  we  can  conclude  from  the  a  priori  error  estimates 
in  Theorem  7  existence  and  uniqueness  of  discrete  solutions. 

Corollary  10  For  a  >  cto  and  (3  >  /30,  the  DG  method  (3)-(4)  admits  a  unique 
solution,  provided  that  h  <  h0. 

Proof  We  need  only  establish  that  if  J  =  0,  the  only  solution  is  u/,  =  0  and 
Ph  =  0.  But  if  3  —  0,  then  u  =  0,  p  =  0  and  the  estimate  of  Theorem  7  implies 
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1  (u h,Ph)  Ilk  <  0  for  h  <  h0.  Since  |||  (•,  •)  \\h  is  a  norm  on  \h  x  Qh,  we  conclude 
that  u/t  =  0  and  ph  =  0.  □ 


The  proof  of  Theorem  7  is  developed  in  section  5  and  section  6.  First,  we  rewrite 
the  DG  method  in  a  non-conforming  fashion,  using  lifting  operators  similar  to  the 
ones  introduced  in  [28],  and  prove  an  inf-sup  condition.  Then,  we  derive  an  abstract 
error  estimate  which  can  be  viewed  as  a  variant  of  Strang’s  lemma.  This  is  done  in 
section  5.  Finally,  in  section  6,  we  make  explicit  the  error  estimates  for  the  principal 
part  of  the  problem  and  then  for  the  T2-norm  of  the  error  by  a  duality  approach 
where  we  need  f3  to  be  large  enough. 


5  Abstract  error  estimates 


In  this  section,  we  prove  an  abstract  error  estimate  for  our  DG  method.  The  key 
ingredient  to  obtain  this  estimate  is  an  inf-sup  condition  which  we  prove  in  Propo¬ 
sition  17,  following  an  argument  often  used  in  the  analysis  of  stabilized  finite  ele¬ 
ments  in  Stokes  flow  (see,  e.g.,  the  survey  article  [43]  and  the  references  therein). 


5.1  A  global  bilinear  form 


For  the  purpose  of  our  analysis,  we  replace  the  integrals  over  interelement  bound¬ 
aries  by  volume  integrals  given  in  terms  of  the  lifting  operators  C  :  L2(£)3  — >-  V/, 


and  M  :  Ll2(£i)  — >  Qh  defined  by 

/n£(v).wrfx  =  J 

f  v  •  {{w  J  ds 

£ 

Vw  G  \h 

J  M(v)qdx.  =  J 

f  vlq}ds 

£z 

Vq  e  Qh- 

We  also  need  the  lifting  operator  ff  : 

L2(£t)  -+  V, 

defined  by 

J  Af(v)  ■  w  dx  =  j 

[  v  [w]jv  ds 

£z 

Vw  <E  Vft. 

Consider  the  forms  ai;ft(-,  •)  and  bylh(-.  •)  given  by 

£tiift(u,  v)  :=  /  VftXu-VftXvdx 

Jn 

-  ja  [C([U]T)  ■  (Vft  x  v)  +  £([v]r)  -  (V* 

friift  (v,p)  :=k2  f  pVh-vdx-  k2  f  -M([v]at)  p  dx. 
Jn  Jn 


13 


Again,  by  integration  by  parts,  we  have 

&tift(v,p)  =  ~k2  f  v  •  V/jpdx  +  k2  f  C(\p]N)  ■  vdx,  V(v,p)  e  V*  x  Qh. 
Jn  Jn 

For  discrete  test  and  trial  functions,  the  forms  •)  and  •)  coincide  with 

a(-,  •)  and  b(-,  •).  However,  this  is  no  longer  true  for  continuous  functions,  due  to  the 
discrete  nature  of  the  lifting  operators.  Nevertheless,  we  carry  out  our  analysis  in 
terms  of  the  forms  «ijft(-.  •)  and  •)  since  they  have  more  favorable  continuity 
and  coercivity  properties  and  take  into  account  the  inconsistency  of  the  forms  by  a 
variant  of  Strang’s  lemma. 

Introducing  the  global  form  Bh(u,  p:  v,  q)  defined  by 

Bh( u,  PI  v,  q)  :=aiift(u,  v)  +  c(u,  v)  +  d(u,  v)  -  k2{ u,  v) 

+  6iift(v,p)  -  bm(u,q)  +  e(p,q), 

we  can  rewrite  the  DG  method  (9)-(10)  in  the  following  compact  form: 

Find  (u h,ph)  G  Yh  x  Qh  such  that 

Bh(uhiPh,v,q)  =  F(v),  (14) 


for  all  (v,q)  G  \h  x  Qh. 


5.2  Stability  of  the  lifting  operators 


The  following  standard  inverse  inequalities  (see,  e.g.,  [38])  will  be  useful  in  the  rest 
of  the  paper. 

Lemma  11  For  polynomials  r  E  Vl(K),  we  have 

1 

||rllo,arr  <  CinvhK2  ||r||o,K  \r\i,K  <  Cmv^^lkllo.rc! 

with  a  constant  C-im  >  0  only  depending  on  the  shape  regularity  constant  n  and 
the  polynomial  degree  t 

We  start  by  establishing  stability  estimates  for  the  lifting  operators. 

Proposition  12  Let  C  and  M.  be  the  lifting  operators  defined  above.  We  have  that, 
for  all  v  E  V(h)  and  for  all  q  E  Q(h), 

ll^dMMIkn  <  Clift ||h  5[v]t||o,£  ||^(Miv)||o,o  <  Clift ||h  5[?]fv||o,f 

||M(!vk)||o,n  <  Ciift 1 1 H  5Iv]iv||o,£x  ll-^,/"([v]iv)||o,n  <  Ciift ||h  3 [v]jv||o,fx 

with  a  constant  C\\^  only  depending  on  n  and  i. 
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Proof  We  prove  the  first  estimate.  Given  v  =  wh  +  w  e  V (h),  observe  that 
[v]t  =  [w/,]'/'  on  £.  By  the  definition  of  the  operator  C  and  the  Cauchy-Schwarz 
inequaiity,  we  have 


I|£(Iv]t)||o,o  =  sup 


<  sup 


In  £(Mt)  -zdx  _  /£[v]T  •  M  ds 

it  ii  —  o  UU 

||z||o,n  zeVft  z  0,0 


h  2  [v]T 

0,8 

hHzlllo,£ 

z  0,0 

Then,  by  using  the  definitions  of  {{ ■ }}  and  h,  and  the  first  inverse  inequality  in 
Lemma  1 1,  we  obtain 


**MII 


2 

0,£ 


-  C  E  ^llzllo,aiv 

KeTh 


<C  £  \\z\\lK 

KeTh 


C\\z 


2 

o.n- 


This  proves  the  first  estimate.  The  other  estimates  are  obtained  similarly.  □ 


5.3  Continuity 

We  can  state  the  following  continuity  properties. 

Proposition  13  There  exists  a  positive  constant  C  only  depending  on  n  and  i  such 
that,  for  cdl  (u,p),  (v,q)  €  V(h)  x  Q(h ), 

|«iift(u,  v)|  <  Gil  (u,  0)  |||h If  (v,  0)  I/,  |c(u,  v)|  <  G|||  (u,  0)  \\\h\\  (v,0)  §A 

|^iift(u,p)|  <  G|||  (u,  0)  I A I  (0,p)  1  h  |d(u,v)|  <  G|||  (u,  0)  \\h\l  (v,  0)  j\h 

|e(p,g)|  <  Cl  (0 ,p)  I,,  1  (0 ,q)  |||ft. 

Consequently,  for  cdl  (u,p),  (v,q)  6  V(/i)  x  Q(h), 

Bh(u,p-,v,q)  <  Gcontl  (u ,p)  |||h  1  (y,q)  |||fc, 

for  a  continuity  constant  GCOnt  >  0  onh  depending  on  n  and  i. 

Proof  Using  the  first  estimate  in  Proposition  12,  we  have 

|«iift(u,  v)|  <||Vft  x  u||0,n  ||Vft  x  v||0,o  +  ||VA  x  v||0,n  ||^(Hr)||o,n 
+  ||Vfc  x  u||0,n  ||£([v]r)||o,n  <  C|||  (u,  0)  |||A|||  (v,0)  ||A. 

Then,  owing  to  the  second  estimate  in  Proposition  12, 

|^iift(u,p)|  <  k2\\Vhp- £([p]jv)||o,n  ||u||0)n  <  C|||  (u,0)  |h|||  (0,p)  |||A. 

The  estimates  for  c(-,  •),  df,  •)  and  e(-,  •)  are  straightforward.  □ 
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5.4  Inf-sup  condition 


We  show  the  stability  of  the  form  B^  in  the  following  two  steps:  we  start  by  prov¬ 
ing,  in  Lemma  14,  a  Garding  inequality  for  the  form  Bf,  in  terms  of  the  seminorm 
|  (-,  •)  \h  and  then,  in  Lemma  16,  a  stability  estimate  for  Bh(u. p:  —^hPi  — p);  by 
combining  these  results,  we  obtain  the  inf-sup  condition  in  Proposition  17. 

Lemma  14  There  exists  a  positive  constant  C  independent  of  h  and  k  such  that, 
for  all  (u,p)  in  V/t  x  Qh, 

Bh{n,p]  u ,p)  >  C |  (u,p)  1 1  -  A;2||u||^, 

provided  that  a  >  Cfft,  where  Cmt  is  the  constant  in  the  estimates  of  Proposi¬ 
tion  12. 

Proof  First,  we  prove  the  following  coercivity  property:  for  all  u  G  V h, 

«iift(u,  u)  +  c(u,  u)  >  C(||Vft  x  u||o;0  +  a  ||li_*[u]T||y .  (15) 

Using  the  arithmetic  geometric  mean  inequality  \ab\  <  ^a2  +  f^b2,  and  the  first 
bound  in  Proposition  12,  we  have,  for  any  5  >  0, 

aiift(u,  u)  +  c(u,  u) 

=  II x  u||o5n  —  2  /  VaX  u  •  £([u]T)  dx  +  a  ||h  i[u]r||^ 

J  fi 

>  (l  -  £)  IIV*  X  u||S,n  -  J  ||£(Iu]r)||S,n  +  a  ||h^[u]T||§,£ 

X  (  l  -  fj  ||V*  x  u||,]  ,,  +  a(l  -  U1)  ||h-i[u]T||^. 

Owing  to  the  assumption  a  >  C'2ih,  we  can  take  C2ift  <  6  <  a  and  obtain  (15).  Now, 
since  Bh( u,  p;  u,  p)  :=  ai;ft(u,  u)  +c(u,  u)  +d(u,  u)  —  /c2||u||o,n+e(p,  p),  the  result 
immediately  follows  from  the  coercivity  property  (15)  and  from  the  definition  of 
the  seminorm  |  (•,  •)  f.  □ 

Remark  15  If  we  replace  the  form  by  its  nonsymmetric  variant  derived  from 
(11),  Lemma  14  holds  true  for  any  a  >  0.  For  the  symmetric  method  in  (9)-(10), 
we  will  assume  throughout  the  text  that  a  satisfies  a  >  C2m  (the  constant  q0  in 
Theorem  7  is  actually  C2i[t). 

Lemma  16  Let  (u,  p)  G  V/,,  x  Qh.  Then  there  exist  positive  constants  C\,  Ci  and 
C;>,  independent  ofh  and  k  such  that 

1  (VaP,p)  Ih  <  Ci  III  (u,p)  II h 

Bh(u,p] -Vhp,  - p )  >  C2k2\\S7 hp\\\n  -  C3\  (u ,p)  \\  -  C3A;2||u||o)n. 
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Proof  Let  us  first  prove  the  continuity  property.  From  Lemma  11,  the  following 
bounds  hold:  ||hi [V^Jiv llo,^  <  C'HVhpHo.n  and  Iloilo,*:  <  Ch^\\ Vp||0,A-. 
Then,  from  the  definition  of  |||  (•,  •)  \\h,  we  obtain 

III  (Vftp,p)  1 1  =  2  k2\\Vhp\\ln  +  2a\\h-^[VhpjT\\l£  +  k2ft  Ilh^V^Ho^ 

+  k2fi  Yl  4l|Ap||^  +  fc27||h^b]Ar||2, 

KeTh 

<Ck2\\Vhp\\l)Q  +  2a\\h-^VhpjT\\l£  +  k2j\\h-^\pjN\\l£ 

<  Ci  I  (u ,p)  \\l 


for  any  u  e  Vh,  with  C\  =  Ci(Ci„v,  ft)- 


In  order  to  prove  the  bound  for  £>/t(u, p:  —  Vpp-  —p),  we  estimate  separately  the 
bilinear  forms  that  are  involved.  We  consider  first  oi;ft(u,  —  V/jf).  From  Proposi¬ 
tion  12: 


«iift(u, -Vhp)  =  j  £([ VhpJT)  'VftXudx 
nci 

^  1 1 V*  —  o  II  V  /  _ 

l0,£ 


>  - 


|h  >[Vfcp]r| 


1, 


2  -2IIV.XU"2 


o,o- 


For  eft,  •)  and  eft,  •)  we  have 


c(u,-Vhp)  = -a  J^  h  1[u]T  •  [V/ip]t  ds 

Ot  ||  1  IT  -r,  |  |  Q  Ot  1  ,  1  n-  -r,  |  |  Q 

>  -^llh  5Ht|Io,£  -  ^||h  HVhp}T\\0,£ 
e(p,  ~p)  =  -k27\KHp}N\\2o,£  ~  «llh_^[Vftp]T||o1g- 

Let  us  consider  now  6i;ft ( - ,  •)•  We  have  again  from  Proposition  12  with  arithmetic 
geometric  mean  inequalities 


fib- 


iift 


(- Vhp,p )  =k2\\Vhp\\ln  -  k2  j^C(\p\N)  ■  Vhpds 


h-2  h-2  (~i2 

>yl|V*-^l|h¬ 


Ijv||o,£ 


and 


-fiiift(u,  -p)  =-  k2  [  u  •  V/jpdx  +  k2  [  £(|p]at)  •  u  ds 
Jn  Jn 


>--k2 1 
“  2  1 


u 


2 

0,0 


h-2  l2 


m 


2 

o,s- 
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For  df,  •),  we  have 

d( u,  VhP)  =  -  k2/3  JV(h[u]jv)  •  Vhp  dx  -  k2/3  fi2K  V  ■  u  A p  ch 

>  —  4  k2(32  C(jft  ||h^  [u] 


IAllo,£ 


KeTh 
12 


YgllVfePlIi.n 

4*2C2/32  AyV-u||2K-^||V^||2n, 


where  we  used  the  estimates  in  Proposition  12,  the  inverse  estimate  ||Ap||0  ^  < 
ChJcWVpWojc  an(i  arithmetic  geometric  mean  inequalities  with  suitable  weights. 
Finally,  we  note  that 

k2(u,Vhp)  >  -4k2\\u\\ln  -  ^||Vftp||gin. 

Adding  together  all  the  contributions  from  the  bilinear  forms  we  obtain  the  result 
(with  our  choice  of  the  weights,  C2  =  |).  □ 


We  are  now  ready  to  prove  the  following  inf-sup  condition. 

Proposition  17  There  are  positive  constants  C\ ,  C2  and  C3  independent  of  h  and 
k  such  that,  for  any  (v,  q)  G  Vft  x  Qh,  there  is  (w,  s )  G  Vft  x  Qh  such  that 

1  (w,  s)  If/,  <  Ci  1  (v,g)  ||ft 

£ft(v,  g;  w,  s)  >  C2(|  (v,  q)  \2h  +  k2\\Whq\\ln)  -  C3/c2||v||2;0. 

Proof  Set  (w,  s)  =  5{v,q)  —  (Vft<?,g),  combine  Lemma  14  and  Lemma 
choose  5  large  enough. 


(16) 

16  and 

□ 


5.5  A  variant  of  Strang ’s  lemma 


We  prove  the  following  abstract  error  estimate  involving  the  residual 

v,g)  =  Bh(u,p\v,q)  -  T{y), 

which  takes  into  account  the  inconsistency  of  the  formulation  (14). 


Theorem  18  There  is  a  constant  C  independent  of  h  and  k  such  that  the  error 
(u  —  lift,  p  —  Ph)  satisfies 


(u-  lift,  p  pft)  I/,  <C 

+ 


,  Jnf  1  (u-v,p-g)  1ft 
|7lft(u,p;w,  s) 

sup 

(0,0)^(w  ,s)£VhxQh 


(w,s) 


+  k  ||u  —  Uft 


o,n 
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Proof  Fix  (v.  q)  e  V/,  x  Qh.  We  split  the  error  (u  —  u/(.  p  —  ph)  into 
(u-  u h,p-ph)  =  (u-  v,p-  q)  +  (v  -  uft,g  -pft)  =:  {(pu,<pp)  +  (£„,£„). 


We  bound  (£u,  fp),  which  we  may  assume  to  be  nonzero.  By  Proposition  17,  there 
exists  a  nonzero  test  function  (w ,  s)  e  V/t  x  Q/*  satisfying  (16)  with  (v,  <?)  = 
(£u!£p)-  We  obtain 


(C&)  |  =  k2\\U\ln  +  *2||Vh^||S,n  +  I  «»,&) 


<  C'|Bh(^„,^;w,a)|  +  C'A:  ||€«||0,n 

<  C'|Bh(v>«,V’P;w,s)|  +  C|fth(u,p;w,  s)|  +  C  /c2||£Jo,n 

<  cr1!  E  +  Cr1iK'*(u’p;w’s)|2 


+ 


III  (W,  8)  fh 

(w,s)  il  +  Ck^Uln 

1\TZh{u,p;w,s)\2 


<  cs-'j  {(pu,(pP)  I! l  +  cs- 

+  CSI  &,,&)  I\l  +  Ck2\\i 


(w  ,s) 


2 

U  1 1  0,0  5 


for  any  5  >  0,  where  we  used  the  definition  of  the  residual  7 Z^,  the  continuity 
of  Bh ,  arithmetic  geometric  mean  inequalities  and  |||  (w,  s)  h  < 

Hence,  the  parameter  6  can  be  chosen  such  that 


(C4)  il<C\\  I l  +  C 


nh(u,p;w,s) 


ill  (w,  «)  III  1 

Since  /c2||£J^n  <  &2||u  -  uh\\ljn  +  A2|| v?M||gjn,  we  have 


+  ck2uu 


2 

o,n- 


(v„^)ll \l 


+  sup 

(0,0)^(w  ,s)eVhxQh 


nh(u,p-,w,s) 


(w  ,s) 


+  A:2 1 1 u  -  Uhllo- 


The  assertion  now  follows  by  applying  the  triangle  inequality  and  taking  the  infi- 
mum  over  all  (v,  q)  e  \h  x  Qh.  □ 


Remark  19  The  result  of  Theorem  18  holds  true  also  in  the  case  /3  =  0.  The 
positivity  of  (3  and  the  stability  induced  by  the  corresponding  forms  will  be  invoked 
in  the  duality  argument  of  the  next  section. 


6  Error  estimates 


In  this  section,  we  make  explicit  the  abstract  error  estimate  in  Theorem  18. 
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6. 1  Approximation  properties 


First,  we  review  the  approximation  results  for  the  L2 -projection,  for  standard  H1- 
conforming  and  for  curl-conforming  Nedelec  operators. 

Lemma  20  Let  w  G  Hl(K),  t  >  0.  Let  II  be  the  L2 -projection  from  H'{K)  onto 
V((K).  Then  for  m  integer,  0  <  m  <  t,  we  have 

I  tt  I  x-f  7  min{^+l,t}— m  n  n 

\w -Uw\m^K  <  ChK  x  |MI t,K- 


Moreover,  ift>  \, 


II  TT  II  ^  .  niin ,^1  2  ||  || 

\\w  -  Uw\\o,dK  <  ChK  2\\w\\ tjK. 

The  constant  C  only  depends  on  k,  £  and  t. 

Proof  For  natural  numbers  t,  the  first  estimate  follows  from  the  classical  Bramble- 
Hilbert  theory  (see,  e.g.,  [38]);  for  t  =  0,  it  is  a  consequence  of  the  stability  of  the 
L2 -projection.  For  non-integer  t,  it  can  be  obtained  by  interpolation.  The  second 
estimate  follows  from  the  trace  theorem  from  L2(dK)  to  Hl{K),  from  the  first 
estimate,  and  scaling  arguments.  □ 

We  also  need  a  standard  H1  -conforming  approximant,  see  [38],  and  a  Clement 
operator,  as  constructed  in  [44] . 

Lemma  21  The  standard  nodal  H1- conforming  interpolant  11  //i  :  [Hl  hl(C)  n 
— >  [Qh  fl  i?o(f2)]  satisfies 

||u;  -  UHiw\\m!n  <  C  hmm{(+1’t+1}~m\\w\\t+1^  m  =  0, 1, 
for  t  >  The  constant  C  only  depends  on  n,  £  and  t. 

This  lemma  is  proved  in  [38]  for  integer  t,  and  can  be  proved  for  non-integer  t  by 
using  the  arguments  in  [45]. 

For  the  Clement  operator  nCi  :  f/0[  (C)  -»  Qh  Cl  F/(J  (0)],  we  recall  the  following 
result  (see,  e.g.,  [44,  pp.  109-111]). 

Lemma  22  There  exists  a  constant  C  only  depending  on  k  and  £  such  that 
(ll™  -  uciwfi,K  +  hK2\ \w  -  naw\\lK  +  hKl\\w  -  Uaw\\ldK)  < 

keT 


Finally,  we  establish  the  next  lemma  for  curl-conforming  Nedelec  operators. 
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Lemma  23  Let  ncuri  be  the  curl-conforming  Nedelec  operator  ( either  of  the  first 
type  [6]  or  of  the  second  type  [7])  into  V/t  fl  i?o(curl;  Q).  Then,  there  exists  a 
constcintC  =  C(K,£,t)  such  that,  for  any  w  e  ifo(curl;  U)fl.f7*(0)3  with  Vx  w  e 


w  -  ncuriw||curi)n  <  C'/imin{^t}(||w||t;n  +  ||  V  x  w||tjn).  (17) 


Moreover,  ifw  also  belongs  to  H( div°;  0), 


h2[w  -  ncurlw]*||o,£l  <  C  hmm^  (||w||tin  +  ||V  x  w||tin) 

E  /ixl|V-(w-ncuriw)||^y  <  (7/imin{^}(||w||t!n  +  || V  x  w||tin). 

KKeTh  J 


Proof  The  first  part  is  proved  in  [45,  Section  5].  For  the  second  part,  denote  by 
UVh  the  L'2 -projection  onto  V/,  .  By  using  the  triangle  inequality,  the  approximation 
results  in  Lemma  20,  the  first  inverse  estimate  in  Lemma  11,  the  Testability  of  Tty,, 
and  the  first  part  of  this  lemma,  we  obtain 

lM[w  -  ncuriw]]v||o,£X 

<c  E  ^(llw-nvftw||^  +  ||nVft(w-ncurlw)||^) 

KeTh 

<c  Y  h2Kmin{e+1’t}\MlK  +  c  E  ||nVft(w-ncurlw)||^ 

KeTh  KeTh 

<  c/.2“”{'+1',)||w||fiI!  +  c||w  -  nCUIlW||2n 

<  +  ||V  x  w|j,,fi)2. 

To  prove  the  last  estimate,  we  integrate  by  parts  and  obtain 

E  ^||V-(w-ncuriw)||^ 

KeTh 

=  -  E  hK  /  vv  •  (w  -  ncurlw)  •  (w  -  ncurlw)dx 

KeTh  Jk 

+  E  hK  f  V  •  (w  -  ncurlw)(w  -  ncuriw)  •  n Kds  =:  7\  +  T2. 

Ke%  dK 

Let  us  first  consider  the  volume  term  7\.  Using  the  fact  that  V  •  w  =  0,  Cauchy- 
Schwarz  inequalities,  the  second  inverse  estimate  in  Lemma  1 1  and  the  first  part  of 
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this  lemma,  we  get 


Ti  <  (  £  hSc\\ W  •  ncurlw||^)  2  ||w  -  ncuriw||o,n 

KeTh 

<  <?(  £  h\\\v  •  ncurlw||^) 2  ||w  -  ncurlw||0,o 

KeTh 

=  c(  £  4IIV  •  (w  -  ncuriw)||^) 2 ||w  -  ncuriw||0!o 

KeTh 

<  Chmm{tt}  (  £  /&||V  •  (w  -  ncuriw)||^)2  (||w||tjn  +  II V  X  w||tjn). 

Ke% 

Similarly,  we  can  bound  the  term  T2  by 

t2<(  £  ^||v-ncurlw||^)2(  £  /iX||w-ncurlw||2aJ2 

Ke%  KeTh 

<  (  Y.  ft'iliv  •  ncurlw||2  A-)i 

KeTh 

(  £  ^llw-nvftw||^  +  ||nVft(w-ncurlw)||2aJ2 

KeTh 

<Cftmi"W}(  £  ft2K||V.(w-ncuIlw)||S,K)i(||w||t,!]  +  ||Vxw||t,!!), 

KeTh 

where  we  used  the  first  inverse  estimate  in  Lemma  1 1  for  the  term  containing  the 
sum  and  proceeded  as  in  the  estimate  of  ||h2[w  —  IIcuriw];v||o,£z  for  the  second 
term.  This  completes  the  proof  of  the  last  estimate.  □ 


6.2  Error  in  the  principal  part 


We  have  the  following  estimate  of  the  residual  in  Theorem  18. 

Lemma  24  Let  (u,  p)  be  the  exact  solution  and  assume  that  V  x  u  G  HS(Q)3  and 
p  G  Hs+1(Et),  for  s  >  |.  Then,  for  all  (v,  q)  G  V (h)  x  Q(h), 

nh{u,  p;  V,  q)  =  jf  [v]r  'fVxu-  nV/l  (V  x  u)}}  ds 
+  k2  [  Mjv  Ip  -  n QhpJ  ds, 

JSx 

where  nVft  and  l\qh  denote  the  L2 -projections  onto  V/>,  and  Qh,  respectively. 
Moreover,  there  exists  C  independent  of  h  such  that 

\nh(u,p-,v,q)\  <  Chmin{e’s}  1  (v,  0)  \\h  (||V  x  u\\ln  +  ||p||£,0)2. 
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Proof  By  straightforward  calculations  involving  integration  by  parts  and  taking 
into  account  the  definition  of  the  lifting  operators  £  and  A4 .  we  have  that,  for  any 
(v,p)  e  V(h)  x  Q(h), 


lZh(u,p-,v,q)  =||Vxu}-  [v]t  ds  —  x  u  •  £([v]r)  dx 
+  k<2  L  W N  ds  -  k2  fnPM({-v]N)  dx. 

Since 

xu-  £([v]t)  dx  =  nV/l(V  x  u)  •  £([v]r)  dx 
J^P M(fv}N)  dx  =  J^YlQhp  M({-v}N)  dx, 

we  obtain  the  desired  expression  for  72-,  (u,  p;  v,  q ). 


(18) 


For  the  estimates  of  the  residual,  let  us  write  72-,(u,  p;  v,  q)  =:  T\  +  T2,  where 
Ti=  f  [v]T -{Vxu-nyjVx  v)J  ds 

T-2  =  k2  f  [v]jv  Ip  -  n QhP}  ds. 

J  Ex 

By  the  Cauchy-Schwarz  inequality,  the  definition  of  the  norm  |||  (-,  •)  ||| ,  and  the 
second  estimate  of  Lemma  20,  we  obtain  the  following  bound: 


Ti<C\\  (v,  0)  1, 

<C\  (v,  0)  I, 


1 

2 


Similarly, 


t2  <  c  1 

< 

0 

>r~ 

E  h^Wp  ^Qhp\\o,dK 

v  KeTh 

<C\ 

(V.  0)  Ik  ( 

{  E  h^{^\\p\\2s+1^ 

K  KeTh  7 

The  estimate  for  72-,  then  follows. 


□ 


We  are  now  ready  to  prove  the  following  error  estimate. 

Corollary  25  Under  the  assumptions  of  Theorem  7,  there  exists  a  constant  C  in¬ 
dependent  of  the  meshsize  h  such  that 

III  (u  u,, p  p,)  If,  <  C'7imm{£’s}(||u||s,n+||Vxu||S)f5+||p||s+i,n)+C'fc||u-u,||0,o. 
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Proof  Consider  the  abstract  estimate  of  Theorem  18  and  bound  the  infimum  by 
I  (u  —  ncurju,p-  limp)  ll/i,  where  IIcuri  is  the  curl-conforming  Nedelec  operator 
and  11// 1  the  standard  if1 -conforming  interpolant  from  Lemma  21.  Taking  into  ac¬ 
count  the  approximation  properties  in  Lemma  21  and  Lemma  23  and  the  estimate 
of  the  residual  in  Lemma  24,  we  obtain  the  result.  □ 


6.3  Error  in  the  L2-norm 


In  order  to  complete  our  error  analysis,  we  need  to  estimate  the  term  k  ||u  —  U/(||0:<> 
This  is  done  in  the  next  proposition  by  a  duality  approach.  The  main  difficulty  in 
this  argument  is  that  we  can  not  assume  any  smoothness  for  the  scalar  potential  of 
the  dual  solution.  To  overcome  these  difficulties  we  will  have  to  make  the  stabiliza¬ 
tion  constant  (3  large  enough. 

Proposition  26  Let  a  >  |  be  the  regularity  exponent  from  Proposition  2.  Fur¬ 
thermore,  we  assume  that  the  exact  solution  satisfies  the  smoothness  assumptions 
in  (13)  with  s  >  |.  Then  we  have 

*||u  -  Uft||0,n  <CX  /iminfa4} I  (u  _  Uft,p  -  Ph)  ||Jh  +  blU^n 

+  C2  I  (u  -  u h,p-  ph)  ||ft, 

with  Ci  =  Ci(k,  i.  k,  (3 ,  Creg.  a,  s )  and  C2  =  C2(k,  t.  Ceii). 

Consequently,  there  exists  /30  =  /30(k.  £.  Cen)  such  that  for  (3  >  (3q 

k\\u  -  Uft||0,n  <  ^  III  (u-  Uft,p-pft)  ||a  +  C'/imin{£’s}||p||s+i!n, 
provided  that  h  <  h0for  a  constant  h0  =  h0(n,  £.  k,  Creg.  o;,  [3. 7,  a,  s). 

Proof  The  proof  is  given  in  several  steps.  We  start  by  introducing  a  suitable  adjoint 
problem  with  right-hand  side  /c2(u  —  U/J,  denoting  by  (z,  ip)  its  solution,  and  we 
express  k2 ||u  -  Uft||o>n  as  the  sum  of  Bh( u  -  u h,p-  ph ;  z  -  IIcurlz,  ip  -  fl aip) 
plus  residual  terms,  with  IIcuri  the  curl-conforming  Nedelec  operator  into  V/,  fl 
if 0 (curl:  Q)  and  na  the  Clement  operator  into  QhC}Hp(Q)  (step  1).  Then,  estimates 
of  the  residuals  (step  2)  and  of  the  bilinear  forms  in  the  definition  of  Bh  give  the 
result  (steps  3  and  4). 

Step  1:  A  dual  problem.  Let  (z,  ip)  be  the  solution  of  the  dual  problem 

V  x  V  x  z  —  k2z  +  k2Vip  =  k2( u  —  U/t)  in  Q 
V  •  z  =  0  in  fl 

n  OO  <19) 

n  x  z  =  0  on  oil 

ip  =  0  on  <9f 1. 
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By  Proposition  1  and  Proposition  2,  we  have 


IMU,n  +  || V  x  z||ffj{2  <  Cie gk2\\u  -  Uh||o,n  2„ 

||0||i,n  <  Ceil 1 1 u  -  uft||0,n 

for  a  regularity  exponent  a  =  cr(fi)  >  We  may  also  assume  that  a  <  1.  The  main 
problem  in  the  following  arguments  is  that  we  can  not  assume  ip  to  be  smoother  than 
belonging  to  Hp(Q). 

The  solution  (z,  ip)  of  problem  (19)  satisfies 

Bh{ z,  -ip',  v,  q)  -  Tlh( z,  -ip]  v,  q)  =  k2( u  -  uh,  v) 

for  all  (v,  q )  G  Hs{Thf  x  Hs(Th),  s  >  §.  Taking  (v,  q)  =  (u  -  uh,-(p-ph)), 
observing  (13),  we  obtain 

Bh{ z,  -"0;  u  -  uA,  -(p  -  ph))  -  1Zh{z,  -ip]  u  -  u -(p  -  p/,))  =  fc2||u  -  uA||o,n- 
Since  B/,(u-uft,p-pft;zft,V»ft)  =  Rh(u,p]zh,iph)  =  -llh(u,p]z-zh,ip-iph), 

for  all  (z/j,  iph)  6V/,x  and  using  the  skew-symmetry  properties  of  Bh,  we  can 
write 


/^2||u  —  U/jIIo  jj  =  £>/i(u  —  Uh,p  —  ph',  z  —  Zh,ip  —  iph)  +  R\+  R2  +  R3  +  R4,  (21) 

with  residual  terms 

Ri  =  -  f  [u  -  uh]T  •  |V  x  z  -  nVft(V  x  z)}}  ds, 

J  £ 

R2  =  k2  f  [u  -  u hJN{ip  -  II Qhip}  ds , 

J  Ex 

R-S  =  -  I  [z  -  zh]T  -  fVxu-  nVft(V  x  u)}}  ds, 

J  £ 

R4  =  -k2  [  [z  -  zh]NP[p  -  n Qhp}  ds. 

J  Ex 

We  define  zh  =  IIcurlz,  so  that  R3  =  0,  and  iph  =  na ip,  with  IIcurl  the  curl- 
conforming  Nedelec  operator  from  Lemma  23  and  nci  the  standard  Clement  oper¬ 
ator  which  satisfy  the  approximation  property  of  Lemma  22. 

Step  2:  The  residuals.  We  estimate  the  residual  expressions  Ri,  R2  and  R4  in  (21) 
(recall  that  R3  =  0).  Let  us  start  with  R4.  The  Cauchy-Schwarz  inequality  and  the 
approximation  properties  in  Lemma  20  yield 

Ri  <  C\\  (u  -  uh,  0)  |||ft(  Y  MV  x  z  -  nVft(V  x  z)||^)  2 

KeTh 

<  Cha (I  (u  -  uh,  0)  IftllV  x  z||ffin. 


i?i  <  Chak2 1  (u-  u/,,0)  |h||u  -  Uft||o,n- 


25 


Thus,  from  (20), 


The  crucial  term  is  R2.  Again  from  the  Cauchy-Schwarz  inequality,  using  the  sec¬ 
ond  estimate  of  Lemma  20  and  the  stability  estimates  (20)  for  the  dual  solution,  we 
obtain 

R2  <  C(k2f3\\h^lu-  uaIjv-Ho,^) 2  (a2/3_1  ^  h^\\tj)  -IiQh^\\ldK)2 

KeTh 

<  Ckfi-* I  (u  u/j,  0)  lUII^Hqn 

<  Ckp~ 2 1  (u-  u/,,0)  1  1 1 u  Uft||0,o, 

with  C  =  C(k,  t,  Ceii).  Finally,  for  R4  we  have 

j?4  <c|(z-ncurlz,o)b(«;2r1  E 

KeTh 

<  c  I)  (z  -  ncuriZ,  0)  |ft|b||s+i,o. 

Step  3:  The  term  |||  (z  —  IIcuriz,  0)  |||ft.  We  claim  that 

I  (z  -  ncuriz,  0)  I/,  <Chak  ||u  -  Uh||0ln»  (22) 

with  C  =  C(k,  £,  k,  Creg,  ft,  a). 

To  see  (22),  we  first  note  that,  since  the  Nedelec  projection  is  curl-conforming, 
||h_s[z  —  ncurlz]T||o,f  =  0.  Furthermore,  by  (17)  and  (20), 

II V  x  (z  -  ncurlz)||on  +  k2 ||z  -  ncurlz||^n  <  C/i2cr(||z||a;n  +  ||V  x  z|k)0)2 

<  CH2ak%n-nh\\l^. 

From  the  second  and  third  estimates  of  Lemma  23  and  (20),  we  have 

k2j3  ||hi [z  -  IIcurlz]]y Ho^  <  Ch2ak2 ||u  -  u^||q,o 

2  4IIV  •  (z  -  ncuriz)||2^  <  Ch2"k2\\u  -  Uft||2n. 

KeT 

The  proof  of  estimate  (22)  now  follows  from  the  definition  of  ||j  (z  —  IIcuriz,  0)  j|jft. 

Step  4:  The  assertion.  We  are  now  able  to  complete  the  proof  Proposition  26. 
Define  =  z  —  ncuriz  and  =  ip  —  IIclb.  From  (21),  taking  into  account 
that,  due  to  the  conformity  of  the  projectors  TTcuri  and  nCi,  c(u  —  Uft,  £z)  =  0  and 
e(p  ~  Phi  iip)  —  0,  we  have  that 

k2 ||u-  Uft | |o>n  =aiift(u-  Uft,^J  +  d(u-  Uft,^)  -  k2( u-  Uft,£J 

+  ^lift (£,ziP  ~  Ph )  —  i’lift (u  —  Uft,  £0)  +  Ri  +  R-2  +  i?3- 

From  the  continuity  properties  of  Proposition  13,  we  obtain 

k'2\\u~  Uft||o,n  <^1  (u  Uft,  p  Ph)  Ilk  III  «„0)  lift 

+  klift(u  —  Uft,  ^0)  |  +  Ri  +  R2  +  R3, 
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where  we  isolated  the  term  &iift(u  —  rqt.  that  needs  to  be  treated  separately.  By 
the  Cauchy-Schwarz  inequality  we  have 


|*iift(u  —  u  ft,^)|5 


[  V/t  ■  (u  -  u/t)  dx  +  k2  f  A^([u-ufc]jv)^<ix 


<C2(k2p  ^4llv h  ■  (u  -  u^Ho,*  +  k2 p\\M(hlu 

v  KeTh 


with  C  =  C(k).  Now,  using  the  approximation  property  of  Lemma  22  of  the 
Clement  operator,  the  third  estimate  in  Proposition  12,  the  second  estimate  in  (20) 
and  the  definition  of  |||  (u  —  uh)  j^,  we  get 

|6|ift (u  - 

<  C^/T^/c2/?  hK\\Vh  ■  (u-  u/t)H^  +  /c2/3||h^[u-  uft]JV||gigI')||^||?in 

v  KeTh  J 

<  C2^^-1 1  (u-  uft,0)  l!l|u-  uj2^, 

where  C  =  C(k.  t.  Cvn).  Inserting  this,  (22)  and  the  estimates  for  the  residuals 
obtained  in  Step  2  in  (23)  completes  the  proof.  □ 


The  proof  of  Theorem  7  follows  now  from  Corollary  25  and  Proposition  26. 


7  Conclusions 


In  this  paper,  we  have  carried  out  an  error  analysis  for  stabilized  interior  penalty  dis¬ 
continuous  Galerkin  methods  for  the  discretization  of  the  indefinite  time-harmonic 
Maxwell  equations.  We  have  derived  error  estimates  that  are  optimal,  provided  that 
the  stabilization  parameters  are  large  enough  and  the  meshsize  is  small  enough.  A 
numerical  study  of  the  proposed  methods  is  the  subject  of  ongoing  work. 
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